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We investigate the quasi-static mechanical response of soils under cyclic loading using a discrete 
model of randomly generated convex polygons. This response exhibits a sequence of regimes, each 
one characterized by a linear accumulation of plastic deformation with the number of cycles. At the 
grain level, a quasi-periodic ratchet-like behavior is observed at the contacts, which excludes the 
existence of an elastic regime. The study of this slow dynamics allows to explore the role of friction 
in the permanent deformation of unbound granular materials supporting railroads and streets. 



A particularly intriguing phenomenon in driven sys- 
tems is the so-called ratchet effect. There is already an 
extensive body of work on this subject, driven by the 
need to understand biophysical systems such as molecu- 
lar motors [1] and certain mechanical and electrical rec- 
tifiers [2] . The classic ratchet is a mechanical device con- 
sisting of a pawl that engages the sloping teeth of a wheel, 
permitting motion in one direction only. Ratchet-like mo- 
tion have been proposed as a mechanism to explain the 
convective motion and size segregation in vibrated granu- 
lar materials [3] . The understanding of this phenomenon 
is crucial in the investigation on the permanent deforma- 
tion in structures subjected to cyclic loading. In soils, 
this loading can be induced by earthquakes, sea waves, 
road traffic, etc. [4]. 

The classical theory of elasto-plasticity describes the 
cyclic loading response by postulating an elastic region 
in the stress space, which changes during the deformation 
[5]. This elastic region, however, is not easy to identify 
because the onset of the plastic deformation is gradual 
and not sharply defined. A great variety of modifications 
have been proposed in order to provide a more appro- 
priate description which, however, make the theory too 
complicated, and require too many material parameters 
which are difficult to calibrate [6] . 

This research was motivated by experiments of cyclic 
loading tests on unbound granular material used to sup- 
port railroads and streets. A slow deformation is ob- 
served during the service life of these structures, where 
the breakage, corrosion and the friction between the 
grains play an important role [4]. We emphasise some 
recent experiments performed in Darmstadt [7]. These 
experiments show that when the samples consist of very 
wear resistant grains, the long time cyclic loading behav- 
ior is given by a linear accumulation of plastic deforma- 
tion. This surprising result suggests that the grains at- 
tain a periodic irreversible motion at the sliding contacts, 
which could, in principle, be detected using numerical 
simulations. Here we report on the first micromechan- 
ical observation using molecular dynamic simulations. 
Ratcheting motion was detected in the sliding contacts 
on a polygonal packing subjected to quasi-static cyclic 
loading. 



The polygons representing the particles in this model 
are generated by using a simple version of Voronoi tes- 
sellation: First, we set a random point in each cell of a 
regular square lattice of side i, then each polygon is con- 
structed assigning to each point that part of the plane 
that is nearer to it than to any other point. This method 
gives a diversity of areas of polygons following a Gaussian 
distribution with mean value £ 2 and variance of 0.36£ 2 . 
The number of edges of the polygons is distributed be- 
tween 4 and 8 for 98.7% of the polygons, with a mean 
value of 6. 

The interaction between the polygons is modelled as 
follows: when two polygons overlap, two points can be 
defined by the intersection of their edges. The segment 
connecting these two intersection points defines the con- 
tact line. Then, the elastic part of the contact force is 
calculated as f e = k n Ax^fi c + k t Ax\t c , where h c and t c 
denote the normal and tangential unitary vectors with re- 
spect to the contact line, and k n and k t are the stiffnesses 
in the respective directions. Ax e n is the overlapping area 
divided by a characteristic length of the interacting poly- 
gon pair [8]. Ax% defines the elastic tangential displace- 
ment of the contact, that is given by the time integral 
starting at the beginning of the contact 
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where O is the Heaviside step function and v% denotes 
the tangential component of the relative velocity v° at 
the contact [9]. This equation defines a limit of elastic- 
ity in the contact force. When the contact force satis- 
fies / t e = /i/^ the contact slides. This produces an irre- 
versible motion between the polygons which contributes 
to a plastic deformation in the assembly. 

In each contact a viscous force f v — mvv c is included. 
Here m is the relative mass of the polygons in contact 
and v is the coefficient of viscosity. This force reduces 
the acoustic waves produced during the loading. Inertia 
effects were excluded by using a very long time of loading, 
indeed, doubling it affects the strain response by less than 
5%. This correspond to the quasi-static approximation. 
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FIG. 1. (a) Shear stress versus shear strain in the first 40 
cycles, (b) permanent (plastic) strain "/m after TV cycles versus 
the number of cycles, (c) stress against the volume fraction 
in the first 40 cycles, (d) volume fraction $jv after N cycles 
versus number of cycles. 

The parameters of our model have been reduced to a 
minimum set of dimensionless parameters: The period 
of cyclic loading to is chosen 4000 times the character- 
istic period of oscillation t s = \Jk n j p£ 2 , where k n is 
the normal stiffness of the contacts and p the density 
of the grains. The relaxation time t r — \jv is chosen 
10 times t s . The tangential stiffness is kt — 0.33fc„; the 
initial pressure po = 0.001fc„ and the friction coefficient 
\i = 0.25. The normal stiffness is chossen k n = WOMPa. 

The simulations are performed on six different sam- 
ples of 400 polygons. First, the polygons are placed into 
a rectangular box such that they do not overlap each 
other. Next, a gravitational field is applied and the sam- 
ple is allowed to consolidate. An external load is imposed 
by applying a force a±W and o^tl on the horizontal and 
vertical walls, respectively. Here <7i and 02 are the hor- 
izontal and vertical stresses. W and H are the width 
and the height of the sample. When the velocity of the 
polygons vanishes gravity is switched off. A fifth-order 
predictor-corrector algorithm is used to solve the equa- 
tions of motion. 

In order to calculate the strain we select the polygons 
whose centers of mass are less than 10£ from the center 
of the sample. Then, the strain tensor is calculated as 
the displacement gradient tensor averaged over the area 
enclosed by the initial configuration of these polygons. 
From the eigenvalues e\ and €2 of the symmetric part of 
this tensor we obtain the shear strain as 7 = ei — £2- 
The volume fraction is calculated as $ = (V p — Vo)/Vb, 
where V p is the sum of the areas of the polygons, Vq the 



sum of the overlapping area between them, and VJ, the 
area of the rectangular box. The initial volume fraction 
is $ = 84.49 ± 0.05%. 

Initially, the sample is isotropically compressed until 
the pressure po is reached. Then, the sample is subjected 
to vertical load-unload cycles as o\ = po + -^Aer[l + 
cos(7rf/£o)] whereas a 2. = po is kept constant. Part (a) 
of Fig. 1 shows the relation between the shear stress 
q = (cti + C2)/2 and the shear strain 7 in the case of a 
loading amplitude Act = 0.6po- This relation consists of 
open hysteresis loops which narrows as consecutive load- 
unload cycles are applied. This hysteresis produces an 
accumulation of strain with the number of cycles which 
is represented by 7^ in the part (b) of Fig 1. We observe 
that the strain response consists of short time regimes, 
with rapid accumulation of plastic strain, and long time 
ratcheting regimes, with a constant accumulation rate of 
plastic strain of around 2.4 x 10~ 6 per cycle. 

Part (c) of Fig. 1 shows the relation between the shear 
stress and the volume fraction. This consists of asymmet- 
ric compaction-dilation cycles which makes the sample to 
compact during the cyclic loading. This compaction is 
shown in part (d) of the Fig. 1. We observe a slow varia- 
tion of the volume fraction during the ratcheting regime, 
and a rapid compaction during the the transition between 
two ratcheting regimes. The evolution of the volume ra- 
tio seems to be rather sensitive to the initial random 
structure of the polygons. Even so we found that after 
8 x 10 3 cycles the volume fraction still slowly increasing 
in all the samples, without any evidence of a saturation 
level. 
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FIG. 2. Plastic deformation per cycle for different loading 
amplitudes . The calculations are performed on six different 
samples. 
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FIG. 3. (a) Each point represents the normal and tangen- 
tial force at one contact. The dotted line represents the sliding 
condition \ft \ — )J.fn- (b) and (c) Trajectories of the contact 
force of two selected sliding contacts. The dots denote the 
times t = 0, 0.5to, 2to in unit of the period to. (d) Plastic 
deformation £ at the contact shown in (c). 

This extremely slow dynamics in the evolution of the 
granular packing shows an astonishing analogy with the 
behavior of glassy systems [10,11]. Even more surpris- 
ing is that no clastic regime is detected by decreasing 
the amplitude of the loading cycles. Fig. 2 shows the 
accumulation rate of strain A7/A./V for different load- 
ing amplitudes Act. A constant accumulation of strain is 
observed during the cyclic loading, even when the ampli- 
tude is as slow as 10~ 3 times the applied pressure. 

The existence of these ratcheting regimes proves to be 
a consequence of the particular disorder of the distribu- 
tion of contact forces. It is well know that the stress 
applied on the boundary is transmitted through an het- 
erogeneous network of contact forces [3] . Both tangential 
and normal force distributions turn out to be very broad. 
The part (a) of Fig. 3 shows the distribution of contact 
forces in a polygonal packing that has been isotropically 
compressed. The components of the force are limited by 
the condition \ft\ < /J,f n , where is the friction coeffi- 
cient. Most of the contacts satisfy the elastic condition 
|/t| < nfn, but due to the heterogenieties, some contacts 
are able to reach the sliding condition |/ ( | = dur- 
ing the compression. Under arbitrarily small load stress 
cycles these contacts slide, given irreversible motion be- 
tween the grains. Those sliding contacts are the most 
relevant micromechanic rearrangements during the cyclic 
loading. Opening and closure of contacts are quite rare 
events, and the coordination number of the packing keeps 
approximately it initial value 4.43 ± 0.08 in all the simu- 



lations. 

After some loading cycles a quasi-periodic irreversible 
motion is observed at these sliding contacts. Two dif- 
ferent sliding modes are represented in Fig. 3: (b) a 
contact sliding forward in the load phase and backward 
in the the unload one and (c) a contact sliding in the 
load phase and sticking in the unload one. A measure 
of the plastic deformation of the sliding contact is given 
by £ = (Ax£ — Axf)/l, where Ax£ and Ax% are the to- 
tal and the elastic part of tangential displacement at the 
contact, the last one being given by Eq. 1. The part (d) 
of Fig. 3 shows the plastic deformation £ of the sliding 
contact shown in (c). Due to the load-unload asymmetry 
of the contact force loop, a net accumulation of plastic 
deformation is observed in each cycle. This is given by a 
slip-stick mechanism which resembles mechanical ratch- 
ets. 

It is interesting to see the spatial correlation of these 
ratchets. Fig. 4 shows a snapshot of the field of plastic 
displacement per cycle at the contacts inside of the speci- 
men. We see that correlated displacements coexist with a 
strongly inhomogeneous distribution of amplitudes. Lo- 
calized slip zones appear periodically during each ratch- 
eting regime. Some slip zones are destroyed and new 
ones are created during the transition between two ratch- 
eting regimes. Moreover, we notice that these ratchets 
are found as well at the boundaries as in bulk material, 
without the layering effects observed in vibrated granular 
materials [3]. 




FIG. 4. The arrows represents the field u of the plastic 
deformations accumulated at the contacts during one cycle: 
u — 500(£at+i — fjv), where £jv is the plastic displacement 
after N cycles. 
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FIG. 5. Fraction of sliding contacts n a in the long time be- 
havior for different values of Aa/po: 0.6 (dash-dotted line), 



0. 1 (dashed line), 0.01 (dotted line) and 0.001 (solid line) 

Let us next discuss the correlation between the dy- 
namics of the sliding contacts and the evolution of the 
stiffness of the material. The last one is given by the 
slope of the stress strain curve in part (a) of the Fig. 

1. The evolution of the fraction n s — N s /N c of sliding 
contacts with the number of loading cycles is shown in 
Fig. 5. Here N s is the number of sliding contacts and 
N c is the total number of contacts. During each loading 
phase, the number of sliding contacts increases, giving 
rise to a continuous decrease of the stiffness as shown in 
part (a) of Fig. 1. An abrupt reduction in the number 
of sliding contacts is observed at the transition from load 
to unload, producing a discontinuity in the stiffness and 
hence, a plastic deformation. We can see also that some 
contacts reach almost periodically the sliding state even 
for extremely small loading cycles. The ratchet-like be- 
havior of these contacts produces a net displacement of 
the hysteretic stress-strain loop in each cycle, ruling out 
an elastic regime. 

The fact that there is no evidence of an elastic regime 
in the cyclic loading response suggests that the hystere- 
sis of the granular materials can not be described by 
using the classical theory of elasto-plasticity. Actually, 
the deficiencies of this theory have been addressed in the 
experimental investigation of soil deformation [6]. Ten 
years ago, the hypoplastic model was formulated in order 
to mend these deficiencies. This continuum theory has 
been supported on the experimental evidence that any 
load involves plastic deformation [12]. 

The hypoplastic approach of granular materials re- 
quires to introduce new ingredients in the current micro- 
mechanical description. The traditional fabric tensor, 
measuring the distribution of the orientation of the con- 
tacts, can not fulfill this description, because it docs 
not make a distinction between elastic and sliding con- 
tacts [13]. New structure tensors taking into account the 
statistics of the sliding contacts, must be introduced in 



order to give a micro-mechanical basis to the granular 
ratcheting. The identification of these internal variables 
and the determination of their evolution equations and 
their connection with the macroscopic variables would be 
an important key to develop an appropriate continuous 
description of granular soils. 

To conclude, we have performed a grain scale investiga- 
tion of the cyclic loading response of a polygonal packing. 
We have shown the existence of long time regimes with a 
constant accumulation of plastic deformation per cycle, 
due to ratcheting motion at the sliding contacts. This 
phenomenon may have deep implications in the study 
of the permanent deformation of soils subjected to cyclic 
loading. More precisely, it may be necessary to introduce 
internal variables in the constitutive relations, connecting 
the dynamics of the sliding contacts with the evolution 
of the continuous variables during cyclic loading. 

At this time, a comparison of the dynamic simulations 
with realistic situations is limited by the computer time 
of the simulations. Using a computer of 2AGHz we are 
able to simulate only 3 x 10 3 cycles per week. Such 
amount of cycles corresponds only to some hours in real- 
istic cases such as railroads and highways [4]. However, 
further investigations of the role of the stiffness, the ini- 
tial stresses, and the inertial effects on this ratcheting 
behavior are currently feasible. 
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